/*********************************************************************************
 * Software License Agreement (BSD License)
 * Copyright (c) 2012, Carlos D. Correa (info@ngraph.org)
 * All rights reserved.
 * 
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. Neither the name of netzen nor the
 *    names of its contributors may be used to endorse or promote products
 *    derived from this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 *  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *********************************************************************************/
/***************************************************************************
 *
 * Eigenvector Centrality Solver
 * based on IETL iterative solvers and Eigen
 *
 * Copyright (C) 2009 by Carlos D. Correa <correac@cs.ucdavis.edu>
 *
 **************************************************************************/
#define EIGEN 1
#ifdef EIGEN
#include "timing.h"

#include "node.h"
#include "centrality.h"

#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <ietl/interface/ublas.h>
#include <ietl/power.h>
#include <ietl/vectorspace.h>
#include <ietl/rayleigh.h>
#include <ietl/iteration.h>
#include <boost/random.hpp>
#include <boost/limits.hpp>

typedef boost::numeric::ublas::matrix<double> BoostMatrix; 
typedef boost::numeric::ublas::vector<double> BoostVector;
typedef ietl::vectorspace<BoostVector> Vecspace;
typedef boost::lagged_fibonacci607 Gen;

#include <Eigen/Core>
#include <Eigen/SVD>

using namespace Eigen;

inline void getDeriv(float *G, float *Qi, int i, int n) {
  for(int k=0;k<n*n;k++) Qi[k] = 0;

  for(int j=0;j<n;j++) {
    if(G[j*n+i]==0) continue;
    float sum = 0;
    for(int k=0;k<n;k++) {
      sum+=G[j*n+k];
    }
    if(sum>0) {
      for(int k=0;k<n;k++) {
		  Qi[j*n+k] = (G[j*n+k]/(sum+1.0) - G[j*n+k]/sum);
      }
      Qi[j*n+i] = ((G[j*n+i]+1.0)/(sum+1.0) - G[j*n+i]/sum);
    }
  }
}

EigenvectorCentralityEigen::EigenvectorCentralityEigen(Controller *controller):AnalysisTool("Centrality","EigenCentralityEigen",controller) {
    setPropertyBool("Derivatives", false);
    setPropertyStringList("Edge Weight", "Default", "edge");
    setPropertyString("Result", this->getName());
	setPropertyStringList("Edge Filter", "Multiplicity", "edge");  /// Temporary. Figure out a general way to apply filters
}

void EigenvectorCentralityEigen::compute() {
  assert(controller);
  string resultProp = getPropertyString("Result");
  string resultDerivs = resultProp + ".Derivative";
  string resultDerivsT = resultProp + ".Derivative.Transpose";
  string resultDerivsM = resultProp + ".Derivative.Mutual";
	string resultDerivsB = resultProp + ".Derivative.Balance";
	string resultDerivsV = resultProp + ".Derivative.Magnitude";
	string resultDerivsN = resultProp + ".Derivative.Negative";
  Graph *g = controller->getGraph("graph");
  bool derivs = this->getPropertyBool("Derivatives");

  compute(g, derivs);
  controller->addDataProperty("node", resultProp, PROPERTY_FLOAT);
  g->computeNodeMinMaxFloat(resultProp);
  if(derivs) {
    controller->addDataProperty("edge",resultDerivs, PROPERTY_FLOAT);
		controller->addDataProperty("edge",resultDerivsT, PROPERTY_FLOAT);
		controller->addDataProperty("edge",resultDerivsM, PROPERTY_FLOAT);
		controller->addDataProperty("edge",resultDerivsB, PROPERTY_FLOAT);
		controller->addDataProperty("edge",resultDerivsV, PROPERTY_FLOAT);
		controller->addDataProperty("edge",resultDerivsN, PROPERTY_FLOAT);
		g->computeEdgeMinMaxFloat(resultDerivs);
		g->computeEdgeMinMaxFloat(resultDerivsT);
		g->computeEdgeMinMaxFloat(resultDerivsM);
		g->computeEdgeMinMaxFloat(resultDerivsB);
		g->computeEdgeMinMaxFloat(resultDerivsV);
		g->computeEdgeMinMaxFloat(resultDerivsN);
  }
}

int EigenvectorCentralityEigen::compute(Graph *graph, bool needDerivs) {
  int n = graph->numNodes();
	printf("Creating matrix of size %d x %d\n", n, n);
  float *G = new float[n*n];
  float *data = new float[n*n];
  float *characteristic = new float[n*n];
  float *qi = new float[n*n];

  string resultProp = getPropertyString("Result");
  string resultDerivs = resultProp + ".Derivative";
  string property = getPropertyString("Edge Weight");
  string resultDerivsT = resultProp + ".Derivative.Transpose";
  string resultDerivsM = resultProp + ".Derivative.Mutual";
  string resultDerivsB = resultProp + ".Derivative.Balance";
  string resultDerivsV = resultProp + ".Derivative.Magnitude";
  string resultDerivsN = resultProp + ".Derivative.Negative";
  print("Saving result as [%s]\n", resultProp.c_str());
	
	string edgeFilter = getPropertyString("Edge Filter");

  if(listener) listener->setProgress(0);

  for(int k=0;k<n*n;k++) {G[k] = 0; data[k] = 0;}
	
	// Build weighted adjacency matrix
	ElementIterator *it = graph->createEdgeIterator();
	for(it->begin();it->hasMoreElements();) {
		Edge *edge = (Edge*) it->next();
		int i = edge->getNode1()->getIndex();
		int j = edge->getNode0()->getIndex();
		float weight = 1.0;
		weight = graph->getAsFloat(edge,PROPERTY(property));
		bool passFilter = graph->getAsBoolean(edge,PROPERTY(edgeFilter));
		weight = (passFilter > 0)? weight: 0;
		// Asymmetric
		G[j*n + i] += weight; 
		// Symmetric
		//G[j*n + i] += 0.5*weight; 	
		//G[i*n + j] += 0.5*weight;
	}  
	graph->destroyIterator(it);

	// Normalize matrix 
  for(int j=0;j<n;j++) {
    float sum = 0;
    for(int i=0;i<n;i++) {
      sum += G[j*n+i];
    }
    if(sum>0) {
      for(int i=0;i<n;i++) {
				data[j*n+i] = G[j*n+i]/sum;
      }
    }
  }

	// Build characteristic matrix Q = I - A
  for(int j=0;j<n;j++) {
    for(int i=0;i<n;i++) {
      double identity = (i==j)? 1.0:0.0;
      characteristic[j*n+i] = identity - data[j*n+i];
    }
  }


  timestamp t1 = now();

  // Initialize Matrix	
  BoostMatrix A(n,n);
  for(int i=0;i<n;i++) {
    for(int j=0;j<n;j++) {
      A(i,j) = data[j*n+i];
    }
  }
  
	// Solve Eigenvalue problem Qx = 0
  ietl::vectorspace<BoostVector> vec(n);
  boost::lagged_fibonacci607 gen;  
  int max_iter = 200; // orig 200
  double rel_tol = 1.0*std::numeric_limits<float>::epsilon();
  double abs_tol = std::numeric_limits<float>::epsilon();
  VectorXf X(n);
  ietl::basic_iteration<double> iter(max_iter, rel_tol, abs_tol);
  	
  printf("[EIGEN] Computing eigenvector Power method max_iters = %d rel_tol = %g\n", max_iter, rel_tol);
  std::pair<double,BoostVector> result = ietl::power(A, gen, iter, vec); 

  BoostVector error = ietl::new_vector(vec);
  ietl::mult(A,result.second,error);
  error -= result.first*result.second;
  std::cout << "error: " << ietl::two_norm(error) << std::endl;

  printf("[EIGEN] Eigenvalue = %f\n", result.first);
  for(int i=0;i<n;i++) {
    X(i) = result.second[i]; 
  }
	
  printf("[EIGENVECTOR CENTRALITY] Time  %f sec\n", now()- t1);
  timestamp t2 = now();
  

	it = graph->createNodeIterator();
	int i = 0;
	for(it->begin();it->hasMoreElements();) {
		Node *node = (Node*) it->next();
		node->setPropertyFloat(resultProp,X(i));
    node->createImplicitEdges(n, graph->nodes);
		i++;
  }
	graph->destroyIterator(it);

  if(needDerivs) {
    Map<MatrixXf> Q(characteristic, n, n);
    MatrixXf Qinv(n, n);
    MatrixXf U(n,n);
    MatrixXf VT(n,n);
    MatrixXf D(n,n);
    MatrixXf T(n,n);
    VectorXf S(n);
	  
    printf("[EIGEN] PseudoInverse of Q\n");

    for(int i=0;i<n;i++) {
      for(int j=0;j<n;j++) {
		  D(i,j) = 0;
      }
      D(i,i) = fabs(S(i))>1e-6? 1/S(i): 0;
    }
	MatrixXf I(n,n);
	I.setIdentity(n,n);
	  
	// Compute pseudo-inverse of Q using SVD
	Q.svd().solve(I, &Qinv);

	// For each node, compute derivative w.r.t each other node
	for(int i=0;i<n;i++) {
		if(listener) listener->setProgress((float) i/(float) n);
		printf("[EIGEN] Computing derivatives %d(%d)\r", i, n);
		 
	// Compute the derivative of the adjacency matrix
		getDeriv(G, qi, i, n);
		Map<MatrixXf> Qi(qi, n, n);	 
      
	// Compute Q^{-1} Q_i X
    // Where X is the vector of centralities and Q_i is the partial derivative of Q w.r.t degree parameter of node i
		VectorXf QiX(n);
		VectorXf Di(n);
	  QiX = Qi*(X*(-1.0));
	  Di = Qinv * (QiX*(-1.0));
		 

	// Store derivatives at edges and implicit edges
		for(int j=0;j<n;j++) {
		  graph->getImplicitEdge(i,j)->setPropertySigned(resultDerivs, Di(j));
		  if(graph->hasEdge(i,j)) graph->getEdge(i,j)->setPropertySigned(resultDerivs, Di(j));

		  graph->getImplicitEdge(j,i)->setPropertySigned(resultDerivsT, Di(j));
		  if(graph->hasEdge(j,i)) graph->getEdge(j,i)->setPropertySigned(resultDerivsT, Di(j));

		  graph->getImplicitEdge(i,j)->setPropertySigned(resultDerivsV, fabs(Di(j)));
		  if(graph->hasEdge(i,j)) graph->getEdge(i,j)->setPropertySigned(resultDerivsV, fabs(Di(j)));
	
		  graph->getImplicitEdge(i,j)->setPropertySigned(resultDerivsN, -Di(j));
		  if(graph->hasEdge(i,j)) graph->getEdge(i,j)->setPropertySigned(resultDerivsN, -Di(j));
		  
	  }
	 }
	  printf("[EIGEN] Computing derivatives DONE\n");
      
	  timestamp t3 = now();

	// Store additional functions of derivatives at edges and implicit edges
	// Such as the mutual derivative and derivative balance 
	  
	for(int i=0;i<n;i++) {
		for(int j=i;j<n;j++) {
			float di, dj;
			di = graph->getImplicitEdge(i,j)->getPropertyFloat(resultDerivs);
			dj = graph->getImplicitEdge(j,i)->getPropertyFloat(resultDerivs);
			float mutual = fabs(di)>fabs(dj)? di: dj; 
			float balance = fabs(di-dj);

			graph->getImplicitEdge(i,j)->setPropertySigned(resultDerivsM, mutual);
			graph->getImplicitEdge(j,i)->setPropertySigned(resultDerivsM, mutual);

			graph->getImplicitEdge(i,j)->setPropertySigned(resultDerivsB, balance);
			graph->getImplicitEdge(j,i)->setPropertySigned(resultDerivsB, balance);

			if(graph->hasEdge(i,j)) graph->getEdge(i,j)->setPropertySigned(resultDerivsM, mutual);
			if(graph->hasEdge(j,i)) graph->getEdge(j,i)->setPropertySigned(resultDerivsM, mutual);

			if(graph->hasEdge(i,j)) graph->getEdge(i,j)->setPropertySigned(resultDerivsB, balance);
			if(graph->hasEdge(j,i)) graph->getEdge(j,i)->setPropertySigned(resultDerivsB, balance);

		}
	}
	  
	printf("[EIGENVECTOR DERIVS] Time  %f sec\n", t3 - t2);
  }
  
	
  
 
  delete G;
  delete data;
  delete characteristic;
  delete qi;
  if(listener) listener->done();

  return 1;
}



#endif //USINGLAPACK
